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ABSTRACT 

We investigate the evolution of interfaces among phases of the interstellar medium with different 
temperature. It is found that, for some initial conditions, the dynamical effects related to conductive 
fronts are very important even if radiation losses, which tend to decelerate the front propagation, 
are taken into account. We have also explored the consequences of the inclusion of shear and 
bulk viscosity, and we have allowed for saturation of the kinetic effects. Numerical simulations of 
a cloud immersed in a hot medium have been performed; depending on the ratio of conductive 
to dynamical time, the density is increased by a huge factor and the cloud may become optically 
thick. Clouds that are highly compressed are able to stop the evaporation process even if their initial 
size is smaller than the Field length. In addition to the numerical approach, the time dependent 
evolution has been studied also analytically. Simple techniques have been applied to the problem 
in order to study the transition stages to a stationary state. The global properties of the solution 
for static and steady fronts and useful relations among the various physical variables are derived; a 
mechanical analogy is often used to clarify the physics of the results. It is demonstrated that a class 
of soliton-like solutions are admitted by the hydro dynamical equations appropriate to describe the 
conduction/cooling fronts (in the inviscid case) that do not require a heat flux at the boundaries. 
Solitons are shown to result from the exact balance of convective and conductive energy transport 
and we demonstrate that they can exist only as associated with the fast velocity mode (analogous 
to the positive Riemann invariant) of the system. Some astrophysical consequences are indicated 
along with some possible applications to the structure of the Galactic ISM and to extragalactic 
objects. 



1. INTRODUCTION 



The problem of the thermal evaporation of cool gas clouds embedded in a hot gas has 
received particular interest in the last years. Apart from the richness of physical aspects involved 
in the process, a strong astrophysical motivation for the study of this type of systems is certainly 
provided by the supposed ubiquitous presence of a hot phase of the interstellar medium (ISM) in 
the galactic disk and halo, coexisting with cooler phases. It is therefore crucial to understand the 
details of the physics to construct plausible models of the ISM. 

The first systematic approach to the problem of the effects of thermal conduction on a cool 
cloud immersed in a hot medium has been undertaken by Cowie & McKee (1977). They obtained 
analytical solutions for the steady-state evaporation flow from the cloud; in addition, they pointed 
out that the classical diffusion approximation of thermal conduction breaks down when the mean 
free path of electrons becomes larger than the temperature scale (saturation). The relevance of 
this effect is quantified by the local parameter ar, which expresses the ratio between classical and 
saturated conduction. Some points are worth to be mentioned: supposing that the flow is isobaric 
(or, equivalently, that its Mach number M is small) implies that classical conduction holds; secondly, 
if conduction is saturated, the usual hydrodynamic equations are not valid anymore, due to the 
large mean free paths involved. Giuliani (1984) pointed out that their solutions are obtained using 
a piece-wise method which leads to unphysical discontinuities of the derivatives of the variables 
at the transition points; his calculation allows instead for a gradual change from the classical to 
saturated regime. 

In a subsequent paper McKee & Cowie (1977) generalized the steady state problem to 
include the effects of radiative losses, and they demonstrated the existence of a critical length at 
which radiative losses balance the conductive heating. A clear cut result for what concerns the 
steady state of conductive/cooling (CC) flows is given by McKee & Begelman (1990). They recog- 
nized that the key parameter governing the structure of the flow is the so-called Field length. This 
length corresponds to the largest wavelength of a perturbation stabilized by thermal conduction 
against thermal instability (Field 1965). Clouds smaller than the Field length will suffer evapo- 
ration, while clouds larger than the Field length will undergo evaporation if the pressure is lower 
than the saturated vapor pressure (Penston & Brown 1970), otherwise they would undergo conden- 
sation. Remarkably enough, this result is in agreement with the one found more than twenty years 
before by Zel'dovich &: Pikelner (1969). These two studies, therefore, clearly define the necessary 
conditions that different ISM phases must fulfill in order to coexist in a stable steady state. 

Aside from the consideration of radiative losses, a number of physical processes may af- 
fect the evolution of conduction fronts. Viscous stresses have been studied by Draine &: Giuliani 
(1984) who show that they may have dramatic consequences on the solutions, in particular when 
conduction is highly saturated. This can be understood recalling that, if ion-electron temperature 
equipartition holds, electron and ion mean free paths are almost equal (Spitzer 1962). The main 
modifications to the results obtained by Cowie & McKee (1977) are a reduced Mach number of the 
evaporative outflow and value of cjt, and the appearance of a collisionless shock transition. They 
also have allowed, in analogy with saturated conduction, for viscosity saturation when velocity 
gradients become large. The impact of conduction on complex multi-phase systems resembling the 
actual ISM has been discussed by Balbus (1985) and Begelman and McKee (1990). 

As it can be realized from the above discussion, most of the studies focused on the steady- 
state properties of the evaporation/condensation process. A remarkable exception is represented by 
the work of Doroshkevich and Zel'dovich (1981). They have studied the temporal evolution of CC 
fronts considering the effect of radiative losses, and have shown that, if the hot gas is freely cooling, 
a cooling wave propagates into the hot phase. Apart from the issue of the uncertain existence 
of the intermediate asymptotic regime they explore (for an extended discussion, see McKee &; 
Begelman 1990), their thermal wave approach neglects all kind of hydrodynamic motions which 
may arise in the system. One of the aims of this paper is to demonstrate the relevance of such 
hydrodynamical effects which, almost unavoidably, are related to the presence of CC fronts, at least 
at the initial stages of their evolution. The basic principle on which our convinction is based is 
that, if the initial gradient of temperature at the interface between the hot and cold phase is large 
enough, the evaporation time becomes smaller than the cloud characteristic dynamical time t^. 
In this case the pressure gradient, driving the motion, relaxes more slowly than the thermal one. 
This importance of the dynamical aspects has been emphasized also by Kovalenko & Shchekinov 



(1992). Thus, the interface behaves somewhat as an explosion site, where the existing pressure 
excess generates an outward flow (evaporation) and an inward one (implosion). 

In a very elegant paper Elphick, Regev & Shaviv (1992) (but see also Elphick, Regev &: 
Spiegel 1991) suggest a new approach to the study of the dynamics of CC fronts. They apply 
non-linear dynamics techniques, as the Lyapunov functionals, to the investigation of the dynamics 
and interaction of a system of localized structures constituted by the fronts in a thermally bistable 
medium. In spite of some approximations introduced into the problem (small thermal conductivity, 
ideal cooling function), their approach reveals extremely promising in the study of the dynamics 
of CC fronts. In the second part of the present paper we also make use of non-linear dynamics 
arguments to obtain the conditions under which a particular class of hydrodynamical non-linear 
waves, i.e., soliton-like, may describe astrophysical CC fronts. Although this subject can have an 
important impact on our understanding of the structure of the ISM (Adams &: Fatuzzo 1992), it 
has not yet been studied in the context of conductive flows. The present paper will concentrate on 
the one dimensional case, but we feel that higher dimensional studies may discover an even larger 
variety of important physical aspects. 

The paper is organized as follows. In §2 we show the result of numerical simulations where 
the effects of radiative losses, viscosity, and ionization are included along, of course, with thermal 
conduction. Saturation effects are also taken into account both for conduction and viscosity. In §3 
we present some analytical calculations that are intended to represent a slightly different formulation 
of some of the aspects of the steady-state analysis, making use of the dynamical ideas put forward 
by this paper; §4 deals with soliton-like solutions, and in §5 a brief summary is given and some 
possible astrophysical consequences of the results are discussed. 



2. EVOLUTION AND DYNAMICS OF CC FRONTS 



Following the line presented in the Introduction, in this Section we mainly want to demon- 
strate that under the conditions which are thought to exist in the interstellar medium of galaxies, 
some regimes may result in strong non-linear, non-steady behaviour of the system. Therefore, 
rather than explore the entire range of parameters of the relevant physical quantities that may gov- 
ern the mass, momentum and energy exchange among the various phases of interstellar medium, 
we concentrate on the cases that can more appropriately show the kind of dynamical features that 
are the main subject of the present work. Moreover, several other studies (see Introduction) have 
already investigated the static and steady-state approximations and we refer the reader to that 
papers for a general discussion. 

The actual interstellar medium is far from being homogeneous and is characterized by 
the presence of strong gradients in the hydrodynamical quantities at the interfaces among the 
various phases by which it is constituted. Even when a pressure equilibrium is achieved between 
spatially contiguous phases, those gradients may induce bulk motions that eventually can drive the 
system away from the initial equilibrium. The mass and energy exchange taking place in these 
cases is governed, in general, by the combined action of thermal conductivity, radiative energy 
losses, viscosity and external heating provided to the system; we can refer appropriately to it as a 
conductive/cooling (CC) front. 

In principle, a magnetic field may also play a non-secondary role in the evolution of CC 
fronts. However, as Balbus (1986) and Borkowski, Balbus & Fristrom (1990) have demonstrated, 
unless the magnetic field is perpendicular to the front propagation direction, it is not likely to 
appreciably modify the general evolution. On the other hand, if magnetic field is tangled and 
threads both hot and cold phases, the reduction of the conductive transport coefficient is not well 
known; some works indicate that this effect can be indeed a minor one, at least for low densities 
(Rosner & Tucker 1989; Tribble 1989). For these reasons, we will neglect magnetic fields in the 
following; the limits of this assumption will be discussed in § 2.1. 

The basic equations describing mass, momentum and energy conservation can be set in the 
following form: 

| = -V.(pv). (2.1) 
p— = -V-U, (2.2) 



-pT— = -pC{p,T,x) + V ■ (kVT) + (n -p<5,,)^, (2.3) 

P = -pT, (2.4) 
/i 

where D/Dt is the Lagrangian derivative. The gas is supposed to be perfect and with solar abun- 
dances, with temperature T, velocity v and density p; the pressure p is obtained from the equation 
of state (2.4); s is the specific entropy; other symbols have usual meaning. 

The equilibrium cooling function C has the usual form 

pjC{p, T) = n2A(x, T) - nH(p, T), (2.5) 

where n is the total particle number density, x is the fractional ionization of the medium and nA 
and Ti arc the cooling and heating rates per particle. We have calculated C from the rates for 
microscopic processes given by Dalgarno & Mc Cray (1972) and Black (1981) at low temperatures, 
while for the optically thin plasma case we have used the results obtained by Raymond et al. (1976). 

The function n is the thermal conduction coefficient governed by atomic diffTision at low 
temperatures (below 10^ K) and by electronic diffiision at high temperatures. The adopted form 
of K is 

k{x, T) = XK^T^I'^ + (1 - x)KaT^I'^, (2.6) 
where Kg = 6.7 x 10~^ and = 3.5 x 10'^ in cgs units. The stress tensor 11 is 

% = -77% - C<^ijV • V +p% 
The non-diagonal components of 11 are 

^''~\dx,^dxi s'^dxij' 

where rj = {2/3)PrK/R is the viscosity coefficient corresponding to shear motions with Prandtl 
number Pr, whereas ( is the bulk viscosity. 

The system (2.1)-(2.4) can be put in a non-dimensional form if the physical variables are 
appropriately scaled. This allows not only to write the equations in a simpler form, but it provides 
a first insight on the physical scales of the problem. A natural scaling is obtained if the cooling time 

Tc = [RT/{^ — l);unA]o is taken as time unit and lengths are expressed in terms of i'p = (kT/ti^ A)o, 
where the subscript indicates that the quantities in parenthesis must be evahiated at appropriate 
values of the variables. The scale length ip is usually named "Field length" ; this length corresponds 
to the maximum wavelength of linear perturbations stabilized by thermal conduction. An additional 
characteristic time = {pi'^/i'j — 1)kT]o, where ^ is a fiducial size of the system, directly related 
to thermal conductivity, can be introduced. 

An important point is that thermal conduction and viscosity are intrinsically kinetic pro- 
cesses and therefore it is crucial to understand the limits of any hydrodynamical approach dealing 
with such phenomena. This limits are set by the condition that the mean free path of the species 
{i.e., atoms or electrons) providing the energy and momentum transport must be larger than the 
scale length of the system, typically represented by the Field length for the cooler gas. As pointed 
out by several papers (Cowie & McKee 1977; Balbus & McKee 1982; McKcc & Begclman 1990) 
this can be a quite common situation in the interstellar and intergalactic medium. A correct so- 
lution of the problem in this case would require the use of Boltzmann equation: unfortunately, 
this is possible only for simple cases where particular simplifications can be introduced (Draine &; 
Giuliani 1984). A good phenomenological approximation is provided by the expression introduced 
by Giuliani (1984) for the thermal flux q 



(2.7) 



(2.8) 



(2.9) 



where gt, known as the saturation parameter, is 



= (2.10) 

Qsat 

and 

Qsat ~ pvlh, (2-11) 

where Vth is the mean velocity for a Maxwelhan distribution of the particle velocities. A similar 
relation holds also for the saturated viscosity (Drainc & Giuliani 1984). In strict analogy with 
conduction they make the following position for the saturated viscosity coefficient 

(1 + (^v) 



where 



2r? 



In the previous equation the parameter represents the largest positive eigenvalue of the matrix 
eij , Pi is the usual ion isotropic pressure and ip^ , describing the maximum temperature anisotropy, 
can be reasonably taken equal to 2. We adopted the same treatment also for the bulk viscosity 
coefficient. 

Before going on and describe the numerical solutions, we believe it is instructive to have 
an insight of the relevant scales of the problem. 

The characteristic lengths of the problem defined above are related by a very simple ex- 
pression: 

2 



(2.13) 

Equation (2.13) describes a line separating the Tc-Tk plane in two regions (Fig. 1): the upper one 
(t/c > Tc) where cooling processes dominate over conductive ones {i.e., condensation) and the lower 
one (rfc < Tc), dominated by thermal conduction (i.e., evaporation). A common interpretation of 
the previous statement is probably based on the implicit assumption that the dynamical behaviour 
of the system can be neglected and the interface between the two phases can be treated as stationary. 
This is equivalent to impose that both Tc and are much shorter than the dynamical response time 
of the system, which can be approximated as ~ {i/cs)o. This characteristic time has in effect a 
slightly more delicate definition when two gas phases with huge temperature differences are present. 
In this case, the actual value of depends on the details of the geometry of the system, but since 
it is likely that the response of the colder phase to external pressure perturbations is slower, its 
dynamical time at least to provides an upper limit to this quantity. An usual astrophysical situation 
consists of a cold cloud, of size i and internal sound speed Cg, embedded in a hot medium, with 
sound speed Cg. To simplify our considerations, we suppose that the cloud is initially in thermal 
equilibrium, £ = 0. If we substitute £ = CgTd in eq. (2.13), the locus of the points satisfying the 
condition Tc = is 

rk=(^Xrl (2.14) 



This curve instersects the one of eq. (2.13) in the point Tc = Ip/cs and, since Tc = Tk on the curve 
(2.13), at the intersection point is = Tc = (Fig. 1). In order to identify the regions of the 
Tc-Tk plane corresponding to a different rfc/r^ ratio, we substitute £ = CgTa in eq. (2.13), with 
the additional condition Td = Tk- These relations provide the locus of the points corresponding to 

Tk = Td- ^ 

Tk=(^-^] T-\ (2.14a) 



These three curves divide the Tc-Tk plane into six regions with different ratios of the three char- 
acteristic times Tc, Tfc, and r^, as shown in Fig. 1. Therefore, once the temperature T of the hot 
medium has been fixed, four different regimes are possible (if we neglect the patologic point at the 
intersection of the three curves), actually depending on £: this cases are labeled in Fig. 1 as A, B, 
C and D. 

Let us consider initially the subregions of this plane for which the relation between the 

characteristic times arc described by strong inequalities (<C or represented by the dashed 
subregions of the Tc-Tfc plane. For a cloud with parameters corresponding to one of these subregions, 
the initial dynamics could be described in terms of the possible evolutionary modes (i.e., isochoric, 
isobaric, isothermal). However, this description is not complete because, as discussed above, three 
different times Tc, , govern the behaviour of the system, and yet it is possible to have different 
intermediate asymptotic regimes sharing a similar initial behaviour. Thus, in order to classify 
different regimes of the evolution of a cloud embedded in a hot intercloud gas it is necessary to take 
into account both the initial behaviour and intermediate asymptotics. Wc will assume, therefore, 
that during the initial stages of the evolution the relation between the largest characteristic times 
is not changed by the different dependence of those times on the dynamical variables. This holds 
for the sufficiently deep parts of the dashed regions of Tc-Tfc plane. The main features of the cloud 
evolution (initially in thermal equilibrium, as mentioned above) can be summarized as follows: 

Region 1 {rk Tc Td)'- initial isochorical heating with subsequent cooling of the heated gas at 
Tc < t <^ Td — no evaporation; 

Region 2 (r^ <C r^; <C Tc): isochorical heating with subsequent expansion of the heated gas at 
Td < t <^ Tc — evaporation; 

Region 3 {Td ^ Tfc <C Tc): isobarical expansion of the slowly heated gas — evaporation; 

Region 4 {Td Tc <^ Tk): isobarical accretion of the cooling intercloud gas close to the interface 
onto the cloud — condensation; 

Region 5 {jc <^ Td <^ Tk): isochorical cooling of the intercloud gas with subsequent slow accretion 

aX Td <t <^Td — condensation; 

Region 6 {jc ^ <C Td)'- isochorical cooling of intercloud gas with subsequent slow accretion — 
condensation. 

A straight consequence is that only the initial evolution of a cloud with parameters lying 
in extreme regions of the Tc-Tfe plane (dashed in Fig. 1) can be considered as static or steady. The 
hydrodynamical motions which develop in the intermediate stages require a dynamical description. 
This statement is even stronger if the cloud initial parameters are located well into the undashed 
parts of Fig. 1, where it is impossible to separate the initial stages of the evolution from the 
intermediate asymptotics. 

The main aim of this paper is to investigate these dynamical effects occurring at the 
interface among different phases of the ISM. 



2.1 Numerical results 

In the following we present the results concerning the numerical solution of eqs. (2.1)-(2.4). 

A number of simplifying assumptions have been made to solve the problem and we want to state 
them clearly to allow the reader to evaluate the limitations of our results. 

We assume a) a spherical cloud immersed in a hot, homogeneous surrounding medium; b) 

initial pressure equilibrium riT = nT, where n and T are the density and temperature of the hot 
gas, respectively, and barred quantities refer to the cloud gas; c) the cloud is supposed to be initially 
in thermal equilibrium (£ = 0) in an isobarically stable state; the hot medium is instead allowed to 
cool down. In practice, given the ratio of the characteristic relevant times for the cases of interest 
here, the temperature of the hot medium is almost constant during the run; d) the ionization 
fraction is calculated at every time step by equating the coUisional ionization and recombination 
rates, e) magnetic fields and non-equilibrium effects in the cooling function are neglected. 

The main difficulty when solving the system (2.1)-(2.4) consists in the large difference 
among the time scales of the various processes. This is a typical situation encountered in stiff 



problems, for which appropriate codes must be developed. For our purposes we have adapted 
the numerical code described by Kovalenko (1993). This is an implicit conservative code in La- 
grangian variables that conserves mass, energy and, in absence of dissipation, entropy as well; the 
conservation of entropy means that the effects of scheme viscosity are negligible. 

Fig. 2 shows the evolution of the relevant hydrodynamical variables describing the cloud- 
hot gas system as a function of radius. For sake of clarity, the displayed variables are normalized 
with the corresponding cloud ones and the velocity is normalized with the cloud internal sound 
speed c<j. The initial conditions are the following: cloud temperature T = 10^ K, cloud density 
n = 1 cm~^, hot gas temperature T = 3 x 10^; the cloud size is ^ = 3.68 pc. From the pressure 
equilibrium condition, thus, the hot gas density is n = 3.3 x 10~^ cm~^. This set of parameters has 
been chosen because they identify a C-type case according to the notation introduced in Fig. 1. As 
already pointed out, these cases are the most suitable one to elucidate the dynamical effects we are 
interested in; by the way, this particular initial condition is also a realistic one as far as the Galactic 
ISM is concerned. Values of the pressure of this order {p/k = 10^ K cm~^) are estimated in the 
Galactic disk, not only for the thermal component but also for the magnetic field and cosmic ray 
ones (Spitzer 1990). In addition, observational arguments indicate that most of the HI mass in the 
thin disk is contained in clouds of sizes 1-10 pc (Lockman, Hobbs & Shull 1986). The Field length 
for such a system is rather large: ip = 1.249 x 10^ pc ^ i, as expected for a dynamical evaporation 
regime. It is also useful to compare the three relevant time scales for this numerical experiment: 
Tc = 3 X 10'' yr, = 8.6 x 10~^Tc, Td = O.OItc. Under these conditions the electron mean free-path 
is Ae = 11.9 pc and therefore thermal conduction is saturated; the global saturation parameter 
(McKee &: Cowie 1977) is do ~ IX^jl = 3.2. The saturation of thermal conduction is included 
in the hydrodynamical code in the approximation described by eq. (2.9). As discussed in the 
Introduction, there are reasons to believe that when thermal conduction becomes saturated, this is 
the case also for viscosity. For sake of comparison, however, we present two numerical experiments 
that explore the cases 99 = oo ( "classical" viscosity) and ip = 2 (saturated viscosity) . 

Six different evolutionary stages, relative to the case = 00, are reported in Fig. 2. At 
the beginning the strong temperature gradient present at the interface produces a huge localized 
thermal energy injection. In some sense, the interface behaves as an explosion site and the problem 
has some analogies, for example, with point explosions. However, the explosion point is not fixed in 
space, but it is comoving with the thermal front towards r = 0. As a result of the energy injection, 
two pressure perturbations are generated, propagating in opposite directions, as illustrated by the 
velocity profiles. The wave travelling with negative velocity becomes supersonic already at about 
t = O.lSr^, while the one travelling away from the cloud surface remains subsonic for all the 
subsequent evolution. This difference can be easily understood recalling that Cg/cs ~ {T/Ty/'^ ~ 
17. The shoulder-like feature shown by the temperature profile in the first three frames has the same 
width of the overpressured region between the two pressure discontinuities. This feature is due to 
the density contrast between the left and the right boundary regions: radiative losses are enhanced 
where the density is higher and therefore the whole zone is strongly radiatively cooled. Another 
interesting phenomenon is represented by the compression of the cloud gas (for which T < 10^ K). 
The spherical compression wave travelling towards the cloud center increases the density of the cool 
material and, at the same time, it is decelerated by the pressure gradient building up at the cloud 
center. Therefore a core is formed with characteristic density ~ 100 times the initial one. We will 
discuss this point in more detail below. When the pressure gradient inside the cloud has reached 
a critical value, the implosion is stopped and the velocity changes sign, i.e., the wave is reflected. 
This is evident from the fifth panel of Fig. 2, where also a pressure spike, which will later relax 
driving the reflected wave, is shown. The reflected wave will eventually superpose to the initial 
outward wave, which at the time t = has already reached r ^ Note that the internal pressure 
is one order of magnitude larger than at t = 0. We have also followed the evolution for t > 
which can be summarized as follows. The mass loss in the outflow is reduced to very small values; 
the system settles down in a quasi-steady state, if one neglects the damped small oscillations of 
the interface around its equilibrium position, occurring with a period of the order of At = 0.25rd. 
These oscillations are the remnant of the initial bouncing and they are the result of the "elasticity" 
of the system. Eventually, at t ~ Tc = lOOr^, the hot medium will start to cool considerably and a 
cooling wave will propagate into the ambient medium. 

An important result is that even if the cloud size is smaller than the Field length, the 



cloud will not be completely evaporated in cases like the present one (grouped under type C in 
Fig. 1) because of the strong compression taking place in its internal regions. In fact, the density 
and pressure increase induced by the compression enhances the radiative losses, thus decreasing, 
and eventually reversing, the front propagation velocity. A similar retarding effect due to cooling 
losses can be realized also from an inspection of the isobaric case, as wc will show in § 3.1, and it 
is responsible for the transition to the steady state. However, the inclusion of the dynamics dra- 
matically amplifies its consequences. In other words, the non-linear effects related to the dynamics 
{i.e., compression waves) modify the ratio between conductive heat flux and cooling, ultimately 
changing the Field length. 

It is also instructive to look at the evolution of the local saturation parameter ax at various 
radii. In Fig. 3 the behaviour of itt is shown at three different evolutionary times. Clearly, inside 
the cloud the conduction is always not saturated. In the interface zone, instead, or can be as high 
as ~ 100, but at t = O.GTd it has already decreased to ~ 8 and after that moment the conduction 
rapidly becomes classical. The peak values of ax are reached closely behind the interface, where 
the temperature gradients are larger; however, secondary structures appear (i.e., the peak at larger 
radius for t = 0.6rd and the dip for t = To) that are mainly governed by the density changes 
corresponding to pressure jumps (sec cqs. [2.10] and [2.11]). 

Fig. 4 shows the evolution for the same parameters adopted in Fig. 2, but now allowing 
for saturated viscosity {(p = 2). The conclusions discussed in the previous case still hold qualita- 
tively. The main difference is represented by the larger velocity of the outward wave caused by the 
decreased efficiency of the viscous dissipation of the pressure gradient. Therefore, being saturation 
of viscosity much higher in the hot medium than in the cloud, the gas compression is much lower 
(~ 10). Saturated viscosity, however, does not seem to provide a mechanism capable to improve 
the cloud evaporation efficiency considerably. 

For completeness, we have performed several other runs for different initial conditions. On 
this basis we can state that the numerical experiment described before has proven particularly 
exemplifying, as long as C-type initial states are considered. In other words, maintaining approx- 
imately the same value for the ratio Td/Tc, but changing the temperature of the hot medium, the 
results are qualitatively the same. 

The cloud implosion process we have shown deserves some more detailed analysis given its 
astrophysical relevance. To provide a more intuitive tool to the discussion, we have reported in 
Fig. 5 (Plate 1) the evolutionary 2-D images of the cloud for the same case presented in Fig. 3. 
The use of density isocontours is aimed to point out the implosion effect. From Fig. 5 the creation 
of a central dense core is quite evident. The maximum value of central compression is reached at 
t = 0.7Td, when n(0) = 253 cm~^ and the cloud size is 1 pc; the temperature in the cloud is in the 
range 1.8 x 10'^ < T < 2.5 x 10^ K. Initially the cloud has a mass of 5 Mq, and at t = O.Ttj it has 
decreased to only 0.19 Mq. This corresponds to an average mass flux from the interface of about 
1.7 X 10~^ Mq yr~^. After this moment, however, the rate of mass loss is rapidly quenched, as 
already discussed, and the cloud is able to survive, even if with a much lower residual mass. 

Another consequence of the implosion process may be of some relevance. The initial column 
density of the gas for the case presented is Nh = 1-13 x 10^^ cm~^; at the moment of maximum 
compression Nh is increased by a factor 7.2. The cloud can become optically thick and, if some 
dust is present, molecule formation may take place in the interior shielded from energetic radiation. 

We turn now to the discussion of the possible model limitations arising from the neglection 
of the magnetic field, B, and non-equilibrium cooling. As far as B is concerned, we can estimate 
its dynamical influence in the following way. When B = 0, the front propagation is stopped, as 
seen above, at the radius at which the internal pressure of the cloud has reached a value Pm, 
roughly equal to the pressure excess caused by conductive heating in the interface. The location 
of this point is determined by the specific characteristics of the dynamical evolution and energy 
balance of the front; in the example studied here, ~ 0.2. When a magnetic field is present, a 
simple but meaningful way to extend the previous analysis makes use of the virial theorem. For a 
spherical cloud, and neglecting thermal pressure, this is written as 

Anr^Pm = (l/^yB^ 

Using the flux-freezing condition irr^B = $ = const., we obtain the equilibrium radius for the 



magnetic (and static) case: 

where the subscript indicates initial values of the quantities. Assuming equipartition between 
magnetic and thermal pressure, from Fig. 2 we derive /3o ~ 0.01; therefore, in the magnetic case 
Tm = 0.28. According to this simple estimation, the main results obtained here should not be 
appreciably modified by the inclusion of B, because its effects would become important roughly 
at the same stages as the thermal ones. This is not surprising if the initial condition prescribes 
equipartition. Of course, these arguments cannot substitute a much needed calculation in which 
magnetic fields are included. As already mentioned, the only other possibility to suppress compres- 
sion requires the unlikely situation of a perfectly insulated cloud such that thermal conduction is 
almost completely inhibited. 

Non-equilibrium effects can certainly be important for the prediction of radiative properties 
of the CC front. These effects have been studied for example by Ballet et al. (1986) and Slavin 
(1989). However, looking at Fig. 1 of Slavin (1989), in spite of large differences in the cooling 
functions, the non-equilibrium temperature profiles obtained differ very little from the equilibrium 
ones. This suggests that the dynamics, mainly governed here by temperature gradients, should 
hardly be affected. Thus, we believe that the equilibrium cooling function used provides at least a 
reasonable first approximation for the study tackled here. 

3. STEADY-STATE PROPERTIES 

The solutions discussed in the previous section illustrate the richness of physical effects 
that may descend from the inclusion of the dynamics in the study of CC fronts. Nevertheless, 
the cases wc have described may appear somehow "extreme", in the sense that the cold phase is 
strongly affected and eroded by the interaction with the hot gas. Such violent dynamical effects 
are driven by huge differences in the energy fluxes transported through the system: these models, 
in fact, correspond to situations in which heat input by thermal conduction and radiative losses are 
strongly out of balance. This naturally leads to the onset of convective motions, as, for example, it 
is well known from the theory of stellar interiors. In another formulation, according to the mapping 
of the Tc-Tk plane previously introduced, we could say that for these states the various time scales 
are very different and thus they are situated far from the bisectrix, where = Tc- Although only 
a minority of the real systems are so finely tuned to be located in the narrow region of the plane 
closed to the bisectrix, they provide an interesting sample of dynamical behaviours. In addition, 
such a bisectrix may represent the final stage of the evolution of a strongly dynamical system, after 
a transient period during which the pressure gradient has relaxed. 

When the pressure can be considered almost constant, the stress tensor reduces to 11 = 
const, and we are left with the energy equation (2.3). This allows us to study the transition from 
the last stages of the dynamical evolution into the steady state. In the following we will derive a law 
for the damped motion of the CC front based on simple flux balance principles. The same approach 
will prove to be quite powerful also in order to understand the soliton-like solutions discussed below. 



3.1 Slowering of the front due to cooling 

For n = const motions are subsonic and the energy equation (2.3) reduces to the standard 
heat equation modified by the presence of cooling. To outline the main physical aspects, we study 
the simple planar case, for which eq. (2.3) becomes 
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(7-l)/x"9t 



d_ 
dx 



dT 
dx 



p£{p,T,x); 



for isochorical processes the previous equation can be set in a non-dimensional form as described 
for the system (2.1)-(2.4) 



dT 

'dt 



d_ 

dx 



dT' 

' dx 



-jCip,T,x); 



(3.1) 



to simplify the notation wc have maintained for n/ p the same symbol k. In order to restrict the 
number of assumptions, we just take H = and assume that A does not depend on the ionization, 
which is a good assumption for temperatures > 10* K. To simplify the notation, we replace C{p, T, x) 
with A(T) in cq. (3.1). The problem then is completely defined when two boundary conditions 
(BCs) are assigned along with the initial condition (the problem is of parabolic nature) 

r(-oo,t) = 0; 
T(l,t) = l; 

T{x,0) = x<0. 

The discontinuity in the initial temperature profile will generate a thermal front, whose position 
we define as the point Xf at which T = 0, propagating towards — oo. A numerical solution of 
eq.(3.1) with the cooling function described in §2 is shown in Fig. 6. The effect of the cooling 
on the dynamics of the front is evident, along with the deceleration of the front to Xf = 0, which 
occurs at i ~ 0.15 and x = —0.36. We interpret these results as the tendency of the system 
towards the state in which conductive flux and volume energy losses are balanced, thus leading to a 
truly stationary state. In order to substantiate this statement we have performed some additional 
analytical calculations. 

From a dimensional analysis of eq. (3.1) we obtain an approximate expression for Xf when 
radiative losses are neglected: 

OT Xf (X t^/^, having taken k cx T^. Note that the usual requirement for a fixed initial energy 
input (Zel'dovich & Raizer 1965) here has been released. It is well known that this problem, with 
different BCs, admits an exact self-similar solution given by the differential equation 



where / and ^ are appropriate self-similar quantities (Zel'dovich &; Raizer 1965). As shown in §2, 
the inclusion of cooling processes introduces different scaling laws. An expression for the shape of 
the front can be obtained from eq. (3.1), looking for solutions of the form T = T{x — vt), thus 
reducing the original partial differential equation to an ODE. The solution, valid for £ = 0, is 

Tlx) ^ f3v\xf -xl^/f^. (3.4) 

To study the transition to a steady state we use flux balance arguments. The energy input £k due 
to thermal conduction can now be dissipated by radiative energy losses £c- Hence the front must 
decelerate and eventually be stop, for a certain value implicitly defined by 

l^fe| = ^ = A(r) = |^,|; (3.5) 
integrating eq. (3.5) over a column of gas of unit area we have 

/ — dx= / X{T)dx. (3.6) 
Jo Jq 

The shape of non- linear conductive front (eq.[3.4]) is characterized by an almost constant profile 
and a steep cutoff close to Xf. Therefore, a good approximation, when evaluating the rhs integral, 
is to assume T = const, which gives 

l^cl ~ x}. (3.7) 



If we use eq. (3.4), and the relation Xf = (l/2)xj'^, we can write the Ihs integrand of eq.(3.6) as 

(3.8) 



dT 1 . .1 



[Xf -X) 



It is easy then to show that 



or, equating to eq. (3.7), 



1 *i-2 

Xf ' 



4(1 + /?) 



x} ~ [4(1 + /?)] 



0/(1-3(3) 



(3.9) 



(3.10) 



For electron conductivity, /3 = 5/2 and x*jr ~ 0.36, in excellent agreement with the numerical 
solution shown in Fig. 6. 



3.2 Static CC configurations 

As demonstrated above, the energy flux drive the system towards an equilibrium state 
in which the energy input (due to thermal conductivity) and output (due to radiative losses) are 
equal. In a CC front, energy is transported by conductive = — kVT and convective qv = pev 
flux, where e is the specific thermal energy. In equilibrium must be 



y"(qv + q«)dS = y P>C • dV, (3.11) 



where the integration is performed in the rhs over the volume of the system, and in the Ihs over the 
surface delimitating the volume, in the appropriate frame of reference in which d/dt = 0. We shall 
consider an cnsamble of scalar variables p, depending on x, p = p{x), and velocity v = (w(x), 0, 0), 

X G (— cxD, +oo). 

To illustrate our method we consider first the simplest case of a static medium, v = 0. 
This method is based on the local analysis of the differential equations governing the system, and 
on the global properties of the energy exchange which determine the characteristics of the flow 
at the boundaries {i.e., ±oo). A similar approach has been developed recently by Elphick et al. 
(1991, 1992) who have investigated the interactions of planar CC fronts. We will concentrate here 
on somewhat different aspects of the problem, namely the general characteristics of the system 
for which steady-state solutions of CC fronts do exist. For the particular case v = 0, the system 
(2.1)-(2.4) reduces to the equation 

= ,£(p.r). (3.12) 



with the BCs 



T(+c^) = Tz; 

the medium is assumed to be in thermal equilibrium at x = ±oo, /2(pi, Ti) = C{p2^ T2) = 0. 

Since in the static case p{T) oc 1/T for isobaric and p{T) = const for isochoric case, 
respectively, we can re- write eq. (3.12) in a simpler form 

A4T)f = AT). (3.14) 

We define the non-dimensional variable 9 = T/T2 and the functions n{9) = k{T)/k{T2), C{9) = 
pC{T)/[p2k{T2)T2], and introduce the new variable r]{9) defined by dr]/d9 = k{9). Given the uni- 
vocal mapping rj 4^ 9, eq. (3.12) can be written also as 



When dC/dr] = 0, it must be \d'^£./dri'^\ > 0. It follows that 

da^ 

= 2>C(r?), (3.16) 

with q = drj/dx = at the two points rj = rji and i] = rj2- The existence of a static solution is 

therefore guaranteed if the medium is thermally stable {dC/drj > 0) at a; = ±oo and Cdrj = 0. 

The last condition on the integral is equivalent to eq. (3.11) and expresses the fact that cooling 
losses are exactly balanced by heating. There is an odd number of stationary points {C = 0) 
pertaining to thermally unstable states with dC/dr] < 0; to simplify our analysis we consider the 
case with a single stationary point, r/o, for which, dC/drjo < 0. 

Near the i-th thermally stable point, the solution of eq. (3.15) can be linearized 

= Cr,{r])Sr], (3.17) 

where drj = rj — rji, Cr, = dC/drj = {dC/dO)/ k. Substituting drj oc exp{u}x) into eq. (3.17), we obtain 

u = (3.18) 

The previous relation describes a saddle point. The characteristic behaviour of the phase trajecto- 
ries near such a stationary point in the q-rj plane is plotted in Fig. 7a-b : 

-1/2 

g~±(£^) \ri-Vi\, for > 0, 

and 

g~(|4j/3)'^'|r?-?7i|3/2, for 4 = 0. 

A peculiar case corresponds to C{r)) ^ \r) — r)i\°' with a < 1, so that = +oo. In this case the 
linearized equation (3.17) is not valid, and a correct solution can be obtained integrating directly 
eq.(3.16), which gives (Fig. 7c) 

q ~ ±\r] - r]i\~- 

The thermally unstable point r/o is a center, and as a result, the phase diagram of the system 
described by eq.(3.12)-(3.13) has the form plotted in Fig. 8a, and the eigenfunction of this system, 
T = T{x) (Fig. 8b), corresponds to the separatrix. It is easy to see that the topology of the phase 
space allows for the existence of a non-trivial solution (T / const) with r(— oo) = r(-|-oo) = T2 
(or Ti) , as shown in Fig. 9a-b. A final possibility corresponds to a periodic structure constituted 
by clouds and intercloud gas (Fig. 9c). 

To clarify the differences among these four cases we recur to a mechanical analogy. We can 
re- write eq. (3.15) in the following form: 

= -V.$(r?), (3.19) 

where = d/drj] $(r/) = — /^^ C{rj')dri' can be seen as a "potential". In the framework of this 

analogy, the characteristic size of the region with T ~ r2 is just the "time interval" Ax needed for 
the mechanical system (3.19) to make a complete orbit around the turning point r] = rj2'- 



Jv' li'n) 



where r)' > rji, for example r)' = rjo. Prom eq. (3.16), the integral is equal to 



Ax{v2)- r J^ ; (3.21) 
J^' V2|$(r?)| 

hence, Ax(?72) is determined by the behaviour of $(??) near r) = 772- Three different possibilities can 
be envisaged: 

a) Let C{r]) = Crj{r] — 772) + J^ririiil — ^2)^/2 + with finite values of CnjCririj ■■■ In this case 
$(77) ^ {rj — 772)^ + ..■ and the integral (3.21) has a logarithmic singularity: 

Ax{r]2 — e) ~ In e, for e — > 0. 

This corresponds to the solution shown in Fig. 8a-b. 

b) £^r]{V2) = 0, £^^(772) / 0; $(77) ~ (77 — 772)'^ + ... and the integral (3.21) has a power-like 
discontinuity at 77 = 772: 

Ax(772 - e) ~ e"^/^. 
Qualitatively, the profile of T{x) is the same as in the previous case. 

c) C{r]) ~ —{112 — 11)"^'^ < 1) so that £^(772) = 00; $(77) r-^ Uo{r]2 — ■ The integral 
(3.21) has a finite value, and the T{x) profile is similar to the one plotted in Fig. 9. 

Thus, static CC fronts corresponding to a i) finite cloud surrounded by hot intercloud 

gas, ii) finite intercloud layer separating two infinite cold regions, arc possible. In addition, a 
configuration constituted by iii) a sequence of flaky cloud-intercloud regions is also possible if the 

cooling function £(77) has a singular behaviour near the stationary points r/i and r/2. These are the 
only solutions for which at least one of the two gas components has a finite size. Note that this 
requires Cj^ ^ and ^(77^) = at the stationary points. 

We conclude that the existence and the structure of steady-state CC fronts is determined by 
the global relation (3.11) and by the local properties of the cooling function and of the conductivity: 
eq.(3.11) expresses the balance of the different mechanisms of energy transport through the system, 

whereas C(T) and k{T) determine the behaviour of the temperature near the stationary (saddle) 
points. An important feature is that, when K(Ti) = 0, the thermal conductivity is able to support 
the energy transport near T = T2 only over a finite interval Ax. In Appendix A we clarify the 
relationship between this conclusion and the Field length. 



3.3 Steady-state CC fronts 

The variety of solutions describing steady-state CC fronts with non-zero velocity is, of 
course, richer than the static one. This is due to the fact that, in general, the interaction of hot and 
cold gas is likely to generate evaporation or condensation flows. Convective motions, in addition, 
lead to a symmetry breaking because, even if the temperature gradient is the same as in the static 
case, the upstream heat flux differs from the downstream one. In terms of the phase space {q,r]), 
this means that a convective heat flux changes the topology of the plane, and the velocity can 
be seen as a bifurcation parameter. Starting again from the integral equation (3.11) wc want to 
obtain the relation between the velocity and the global temperature and density distribution for 
the steady-state case. 

Let us consider a steady-state flow corresponding to an evaporation front travelling with 
velocity uq. All the hydrodynamical variables can be written in the form a = a{x + uot). In the 
frame of reference comoving with the evaporation wave the appropriate equations are: 

V{pu) = 0, (3.22) 



puVu + Vp = 0, 



(3.23) 



-^puVT = -pVu + V{kVT) - C{p, T), (3.24) 
7-1 

p = -pT. (3.25) 

Here u is the velocity in the comoving frame reference, V = d/dX, X = x + uqI. We introduce the 
following non-dimensional variables 



i = X/iF^, p = p/p^,T = T/Tu 

P = 7P/(PiCi),C/ = u/ci, A = Ip^/iiclpi), 



(3.26) 



where the subscript 1 refers to the unperturbed state at X = —00, and Ci = R^yTi/ p,. 
The first integrals of motion are the following: 

pU = Uo = J, (3.27) 

pC/2+p = [/2_^1, (3.28) 

^{pv - 1) + i(f/^ - C/o^) = (3.29) 

where J = Uq is the mass flux, v = p~^ is the specific volume, and 

A,(0= A(p(0,T(0)de (3.30) 
•J —00 

Writing this system for ^ = +00, and adding the equation of state 

A(p,T) = (3.31) 

for the gas at ^ = +00, we can obtain the values of p, p, U at ^ = +00; Uq at ^ = —00 is equal to 
the velocity of the evaporation front for given values of pi,pi,Ti at ^ = —00 and r2 at ^ = +00. 
Note that equation (3.31) is not valid in the intermediate range —00 < < +00. 

Let us assume that X{p,T) = has a shape similar to a van dcr Waals-typc equation, 
shown in Fig. 10, usually appropriate for a two-phase ISM model (Field, Goldsmith & Habing, 
1969). We stress, however, that the results obtained in this Section do not depend on the specific 
form of the cooling function adopted but only on the number of its zeros. Thus, 



p 



{cj/j)p for ^ = -00; 
(ci/7)P for ^ = +00, 



where C2 ^ ci for the differences between warm (T2 ~ W^K) and cold (Ti ~ lO^K), or hot 
(T2 ~ lO^K) and warm (Ti ~ lO^K) ISM phases. In non-dimensional form the equation of state 
for the hot gas (^ = -|-oo) is 

P2 = clv^\ (3.32) 
The dynamical integrals of motion (3.27)-(3.29) can be reduced to 

P2 = 1 + 7U^{1-V2), (3.33) 

and 

^ -iP2V2-l) + k^-p2)v2-^U^il-V2) = Q, (3.34) 



where Q = —Xs/Uo, Xs = As(+oo). Equation (3.33) is sometimes referred to as the Mikhelson 
equation (Zel'dovich &; Kompaneets 1955), while equation (3.34) is similar to a Hugoniot adiabatic 
for As = 0. Combining equations (3.33) and (3.34) with the equation of state for the hot gas (3.32), 
we obtain ^ 

^ = l + jU^(l-V2), (3.35) 

and 

-^(^^ - 1) + |(1 - - ^(1 - V,) = Q. (3.36) 

Equation (3.35) has a solution for jUq <C c^^ 

V2 = cl{l + jU^cl + o(7C/o'c2)); (3.37) 
after substitution into eq.(3.27) we obtain 

U2 = clUo{l + -/U^4 + o{-/Uocl)). (3.38) 

The assumed approximation 'jUq ^ c^^ corresponds to a subsonic flow at ^ = +oo. Under such 
conditions eq.(3.36) can be reduced to 



clUo + 7U^4 = -2X,, (3.39) 



7 

whose solution is 



/ ^2 



For As > cq. (3.40) describes an evaporation wave (net heating, Uq > 0), whereas for As < 
it describes a condensation wave (net cooling, Uq < 0). It can be shown (Appendix B) that this 
represents a subsonic flow at ^ = -j-oo in accordance with the assumption made above. 

We now concentrate on the properties of the phase plane (Hayashi 1985 provides an excel- 
lent background on this topic). It has to be pointed out that the flow considered is approximately 
isobaric, as derived from eq. (3.40): 

p = pc^ = l- ^U^c" + oi^U^c") (3.41) 

This leads to the following form of the energy equation (3.24) 

I ^dT dU d dT ^, 

-Uo-7T + -T-- -jtK— = -A(p, T). (3.42) 



7-1 dC 

Combining the mass and momentum integrals, (3.27)-(3.28), with the equation of state (3.25) we 
find 

1 -I- r/2 

T = -U'' + =-jr^U, (3.43) 
which gives an estimation of the temperature valid for subsonic flows 

After substitution into eq. (3.42), we can re- write the same equation as 

f-^-o|^A„.,, ,3.4.) 



In analogy with the static case, we wiU assume that X{p,rj) has three stationary points (Fig. 11a): 
two of them correspond to thermally stable states = — oo and ^ = +00) whereas the third 
= Co) corresponds to a thermally unstable one. For small perturbations around the stationary 
points we have: 

£5r] J UodSr] 
7 - 1 K 

where Sp ~ —p^drj/n for isobaric perturbations. Taking 5r] oc exp{io^) the characteristic equation 
corresponding to (3.45) is 

and the eigenvalues are 



7 — 1 K \ K 




2 7 — 1 K 

Since for a thermally stable gas p^ /nXp — A,, < 0, it is easy to see that the nature of the ther- 
mally stable points corresponds to a saddle point (Fig. lib). At the intermediate point 770, 
instead, p^/nXp — X-q > 0, When condition (3.40) holds, the discriminant is negative because 
C/2 = o((^f/4)^C2^), whereas C/q = o((£f/4)c2"^) and Ap,A^ = o{If/Q. Thus ?7o is an unstable 
spiral point (Fig. lib), and the phase diagram differs qualitatively from the static case. The main 
difference is that evaporation (Fig. 11c) and condensation (Fig. lid) waves are not symmetric 
anymore. Thus, a steady-state CC front can exist only as a switch-on wave between two stable 
phases (cold/hot), and periodic structures like those plotted in Fig. 9c are not realisable. 

To better understand this difference we apply again the mechanical analogy used above. 
Equation (3.44) can be written as follows: 

(3.47) 

In this case the potential $(r/) is an asymmetric function of rj (Fig. 12a-b); therefore, for the 
evaporative case {dij/d^ > 0) a negative "friction force", like the one represented by the second 
term on the Ihs of eq. (3.50), increases the "energy" of the "particle" during its motion from r] = r/i 
to 77 = rf2- Since $(?7i) = 0, we have $(772) > ^{vi) fo^' an evaporation wave, and ^(772) < ^{vi) 
for a condensation wave. This means that it is impossible for this "particle" to reach = r/2 
from rj = rji and turn back, because such transitions lead to a difference between the initial and 
final values of the "potential" at rj = rji. Out of analogy, this is due, of course, to the radiative 
losses occurred in the gas during the phase transitions. The cigcnsolution of eq. (3.47) corresponds 
to the equality between the "potential" barrier and the "work" made by the negative "friction 
force". Indeed, multiplying eq. (3.47) by drj/d^ and integrating over ^(— oo,-|-oo) with the BCs 
{drj/d$,)±oo = 0, we obtain 

''^ 7 Uq f drjX _ J $(772) — $(771), for evaporation; 



37 dr] 



7 — 1 K \d^ J I $(r/i) — $(772), for condensation. 



An artificial way to maintain a steady-state CC front with a structure of a hot (cold) 
gaseous layer surrounded by cold (hot) gas (similar to that of the static CC configuration plotted 
in Fig. 9) is to inject a fixed heat fiux at a certain ^. In such a way, the additional heat fiux 
behind the condensation wave provides heating to the cooled gas and evaporates it, restoring the 
initial state (Fig. 13a). The specific value of the heat fiux necessary to maintain such a CC front is 
determined by the value of drj/d^ at the top of the phase trajectory coming out of the critical point 
77 = 772 (in Fig. 13a it corresponds to the point S). In the case of an evaporation wave, a negative 
heat flux behind the CC front dilutes the excess of thermal energy, provides a cooling mechanism 



to the hot gas and favours the condensation process (Fig. 13b). However, this requires that the 
heat source moves together with the CC front. 



4. EXISTENCE OF SOLITON-LIKE SOLUTIONS 

The previous results demonstrate that a CC front structure similar to a soliton (conden- 
sation - evaporation [CE] wave or evaporation - condensation [EC] wave) is possible if supported 
by an artificial heat source. However, the hydro dynamical integrals contain non-linear terms that 
may cause the same effect. According to eq. (3.43) 

1 -I- TI^ 

AT = -2UAU + —j-^AU. 

The two terms on the rhs produce a perturbation of T in opposite directions. This perturbation 
is small only if \Uo\ -C 1, as in the static case. Thus, we expect that the non-linear term [/^ in 
eq. (3.43) may result in a qualitative change of the phase plane topology that allows soliton-like 
solutions without any artificial heat source at the boundary. 

Let us consider the energy equation (3.24), describing a stationary wave with velocity Uq 

Uo (IT dU d dT ^, ^. 

-P-jy + ZJ7^-jy-Hp,T), (4.1) 



7 - 1 d^ d^ 

and with BCs r(— oo) = r(-|-oo) = Ti. Note that Uo in eq. (4.1) is in general different from the 
Uq determined in the previous Section for evaporation or condensation waves. Using the relations 
among T, p and U determined by the first integrals we have 

p = U^ + l-UoU = eo- UoU, (4.2) 

T=-U^ + ^^U = -U^ + ^U, (4.3) 

Uo Uo 

where eo = 1 + C/q is the specific enthalpy. Equation (4.3) allows us to reduce eq. (4.1) to a form 
containing only T. For this purpose we solve eq. (4.3) to find U = U{T): 




liWl-^Tl ^/±(r), (4.4) 



where it correspond to the "fast" and "slow" velocity mode, respectively. After substitution of eq. 
(4.2) into eq. (4.1), using equation (4.4), we can re-write eq. (4.1) as follows 



where dQ = d^/n, and 



dC 2 Vt-I J dC 



F^ = ±ll-^T 



-1/2 



We further assume that the cooling function has two equilibrium points: Ti = 1, and T2 
defined by X{p,T) = 0. These two points are stationary points on the phase plane T-{dT/d(), that 
are {Ti,0) and (T2,0). The solution corresponding to a solitary (EC or CE) wave is possible if one 
of the points is a center and the other one is a saddle point. The characteristic equation for linear 
perturbations ST oc exp(a;C) for the two stationary points is 



2 Uo fj + l „±\ / dL\ 



where i =1,2 and dL/dT = K{Ti){dX/dT)i. Hence, the eigenvalues are 



UJ 



1,2 



4 
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dTj, 



(4.7) 



The thermally stable point with {dL/dT)i > is a saddle point, whereas the nature of the thermally 
unstable point with {dL/dT)i < depends on the value of (7+I/7 — 1 — F^^ . Wc analyze now 
in detail the case that corresponding to soliton-like solutions which arc the subject, of this Section. 
To be specific, we suppose that the saddle point is identified with T = Ti, and T = r2 > Ti is the 
center. The phase diagram of the system in this case is plotted in Fig. 14a; the relative temperature 
eigenfTinction is shown in Fig. 14b. 

It is easy to see that the point (T2 , 0) can be a center only for the fast mode solution: 



(4.8) 




with the necessary condition — F+ = 0, equivalent to the equation: 



1- li^Ta 



7+1 

7-r 



(4.9) 



whose solution is 



^0 



(7 + 1)^ 
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(7 + 1) 
27 



(4.10) 



This solution determines two values of C/g > if the sound speed at the center point satisfies the 
inequality C2 > 47/(7 + 1)^) valid for T2 > Ti = 1. In this case, there exist two different soliton- 
like solutions corresponding to supersonic and subsonic solitary waves of EC type. The FWHM 
of the soliton can be estimated directly from eq.(4.5). At the top of the soliton dT/dQ = and 
d^T/d(^ < 0, kA 7^ 0, hence 



Ax 



WT/dQ^ 



All quantities should be taken for T = T^, where T^, is the temperature at the soliton top (or bottom 
for a CE soliton); can be estimated at the point where $ = 0, ($ is defined by eq. [3.19]). 

A similar analysis for the case with a saddle point at T = T2 (thermally stable phase) and 
a center point at T = Ti demonstrates the existence of a CE soliton-like solution with a reversed 
shape with respect to the EC case. 

These soliton-like solutions correspond to the separatrix of the phase plane (see Fig. 14a) 
like the hydrodynamical solitons described by the Korteweg-de Vries equation. However, their 
physical nature is quite different. In fact, they are the result of the simultaneous action of four 
processes: thermal conductivity V(/{VT), the cooling C{p,T), convective transport puVT and 
advection pVu. Indeed, the thermal conductivity determines the second-order characteristic equa- 
tion (4.6) and the set of solutions near the stationary points; the radiative cooling determines the 
number of stationary points and the local topology of the phase space in their vicinity. The most 
interesting role is played by the last two processes: their conspiracy is able to restore the symmetry 
of the phase space broken by the hydrodynamical motions for slow isobaric flows (e.g., eq. [3.33]). 
The sum of the hydrodynamical terms in eq. (4.1), taking into account eqs. (4.2) and (4.3), can 
be written as: 

Uo dT dU _Uq h + l eo dU\ dT 
T^'d^ '^^'di ~ T \^ + TTodTj d^- 



Thus, the first term on the rhs represents an asymmetric heat flux and hence destroys the sym- 
metry of the phase space present in the isobaric case {dU / dT = 0) . The second term describes a 
stabihzation or amphfication of the convective heat flux depending on the sign of dU/dT. For the 
fast mode dU /dT < 0, thus the convective heat flux is stabilized and balanced for an appropriate 
value of Uq at a certain T = Ti (or T = T2). This provides the soliton-like shape of the wave. The 
soliton temperature profile, however, is slightly asymmetric as it can be realized by rewriting eq. 
(4.5) in the following form 



It is clear that the upper and lower parts of the (T, q) plane are asymmetric: at the points where 
^ = 0, in fact, T{q) for g > differs from T{q) for q < 0. 



In this paper we have studied the dynamics of thermal conduction fronts in a multi-phase 
medium with parameters suitable to describe the general ISM. The inclusion of radiative losses 
affects substantially both the dynamics and the structure of the conductive/cooling front. However, 
a number of additional physical effects like viscosity, ionization, saturation of the kinetic processes 
have been shown to be able to modify appreciably the overall behaviour of the system, when 
included. The presence of several, often very different, scales in the problem determines regimes in 
which the response of the system is strongly time dependent and non-linear. CC fronts may induce 
relevant dynamical effects if the conductive time is much shorter that the sound crossing time of 
the cloud. This condition is, in practice, equivalent to require that the cloud size is smaller than 
the Field length for the hot gas. 

As an example, we have explored in detail a realistic dynamical regime by means of nu- 
merical simulations. The adopted initial condition consists of a spherical cloud in equilibrium with 
a surrounding hot medium at a pressure p/A; = cm~^ K, with a ratio of intercloud-to-cloud gas 
temperature of T/T = 300. This set-up may fairly well reproduce a typical interstellar Galactic 
environment. The main results are summarized as follows: 

1. For the case (p = co (classical viscosity) a huge pressure gradient is created at the interface 
hot/cool gas which drives an outward and an inward flow. In some sense, the interface 
behaves like an explosion site. The supersonic wave travelling towards the cloud center 
compresses the cloud and at the same time tends to evaporate it. A dense (> 100 times 
the initial density) core is formed; the subsequent increase in pressure is able to stop the 
compression wave and to reflect it back. 

2. We have shown that after the pressure gradient has relaxed, the front is decelerated, and 
the system tends to a steady-state. Using simplifying positions as neglecting the viscosity 
and thermal conduction saturation, the final stages of this process have been investigated 
exploiting flux balance arguments. These findings have been substantiated by their agree- 
ment with exact numerical solutions. 

3. Clouds of size smaller than the Field length can be able eventually to stop the evaporation 
process due to the density and pressure increase following the compression phase; the net 
results is that the enhanced radiative losses decelerate the CC front (non-linear effect). 

4. In the first evolutionary stages, thermal conduction is highly saturated with a value of 
(Tt ~ 100, but it rapidly becomes classical after about one half of the dynamical time. The 
role of the dynamics in decreasing the saturation level of the kinetic effects is an important 
piece of information to add to the results for the steady cases obtained by McKee &; Cowie 
(1977) and Draine k Giuhani (1984) 

5. If saturated viscosity is allowed, (/? = 2, the compression factor is decreased by a factor of 

10, but the main features of the flow remain qualitatively the same. 

The global properties of static and steady-state CC front solutions have been derived from 
an analysis of their behaviour in the phase plane. Different thermal phases individuate stationary 
states whose topological nature is determined by the stability properties of those states. Depending 




5. SUMMARY AND POSSIBLE APPLICATIONS 



on the behaviour of the coohng function at the stationary points, it has been shown that in the 
static case three different configurations are possible: i) finite cloud surrounded by hot intercloud 
gas, ii) finite intercloud layer separating two infinite cold regions, and iii) a periodic sequence 
of cloud-intercloud regions. The last possibility is related to the fact that in this case thermal 
conduction can transport energy only over a finite spatial interval. The addition of a velocity field 
produces a symmetry breaking in the topology of the solutions and, therefore, periodic structures 
are not allowed. 

It has been shown that solitary wave solutions arc admitted by the set of equations describ- 
ing CC fronts, even in the case of no heat flux at the boundaries, thanks to the intrinsic non-linear 
character of these equations. Physically, those soliton-like solutions are the result of the exact bal- 
ance between convective and conductive heat flux, which restores the symmetry breaking occurring, 
in general, when the flow is not isobaric. These solitary waves have the remarkable property that, 
differently from an evaporation front that leaves the gas at a higher temperature behind it, they 
change the temperature of the medium affected by their passage but after this transient the gas 
returns to its initial state. 

Finally, we add here that it is possible to demonstrate the existence of self-similar solutions 
of the equations describing CC fronts in absence of bulk motions. Their properties are currently 
under investigation and they will be discussed elsewhere. The preliminary results show that, con- 
trary to the suggestion of Doroshkevich k, Zel'dovich, cooling media do not necessarily induce 
condensation waves, but they may excite evaporation modes in the system. 

The most severe limitation of our model descends from the spherical geometry used in 
the numerical calculations. The main reason for this choice is, of course, dictated by simplicity 
requirements; real ISM clouds can largely deviate from this idealisation. Some of the findings of 
this paper can be modified by more refined and multi-dimensional calculations; the striking nature 
of the implosion effect, in particular, can be partially an artifact of the geometry. However, we 
believe that the basic physics (compression, density increase, slowering of the front, wave reflection) 
underlying the process, and outlined throughout the paper, must be realized also by more complex 
configurations. In addition, this simplification highly facilitates the comparison with analogous 
classical works in the field, based on the same assumption. 

Some mathematical problems related to the existence of solitons need some more accurate 
analysis. For example, it is not clear to us if the the non-linear terms present in the hydrodynamic 
equations are likely to destroy closed trajectories in the phase space around the center point for 
large amplitudes. In addition, the possibility of closed paths seems to depend on the nature of 
the system, which in turn is determined by its characteristic equation (Guckenheimer & Holmes 
1983). Finally, the stability of solitons for the Korteweg-de Vries equation has been demonstrated 
by Benjamin (1972), but the analogous question for the thermal solitons we are presenting here is 
completely open. 

Although this paper is not particularly dedicated to a specific astrophysical situation, a 
number of consequences worth to be explored can be suggested. 

A result that has important and immediate implications for our understanding of the ISM 
descends from the cloud implosion induced by thermal conduction. Due to the unavoidable density 
increase, the cloud tends to become optically thick, providing a site for molecule formation. A clear 
evidence of molecular clouds embedded in neutral hydrogen, both surrounded by hot, x-ray emitting 
plasma, has been provided by the ROSAT shadowing experiments (Burrows & Mendenhall 1991; 
Snowden et al. 1991) in the direction of the Draco complex. The phenomenon we have studied 
requires cloud properties which are by far typical of the HI clouds populating the Galaxy. Therefore, 
it may take place whenever very hot gas (T 3 x 10^ — 10'' K) is suddenly injected by some process 
the most natural being supernova explosions) in their surroundings. In a couple of papers, Bertoldi 
1989) and Bertoldi & McKee (1990) have discussed an analogous cloud compression mechanism 
based on the action of an ionization-shock induced by the radiation field of newly born stars. 
Our mechanism is rather important in the regions in which the radiation field in not particularly 
enhanced, but some hot gas is present. A prototype of these regions can be the Galactic halo, where 
hot gas is suggested by several different observations (Sembach & Savage 1992; for a review see 
Spitzer 1990) along with neutral HI clouds. Those clouds can be a product of a thermal instability 
occurring in a fountain flow under particular conditions (Houck & Bregman 1990, Ferrara & Einaudi 
1992, Li & Ikeuchi 1992). However, it has been pointed out (Ferrara &: Einaudi 1992) that for the 



prevailing conditions of the hot flow, turbulence is more easily generated than a condensation; in 
this case the clouds should have a different origin. If the predictions of the McKee &: Ostriker 
(1977) model are correct, suitable conditions for the implosion mechanism can easily be found in 
the disk as well. 

Some other aspects of the cloud/hot gas interaction deserve additional study. For example, 
clouds engulfed by a strongly anisotropic hot gas flow can be rocket-accelerated as a whole, by the 
same mechanism providing implosion in the spherically symmetric case. This may have important 
consequences for the observed degree of HI turbulence required by some models to support the 
extended neutral Galactic component (Lockman & Gchman 1991; Ferrara 1993). Also, the inter- 
action of a cloud with the hot surrounding medium may induce turbulent motions inside the cloud. 
These motions are likely to be supersonic, as can be realized from the results of §2.1. Therefore, 
even in presence of a substantial dissipation, the energy flux due to conductivity may be able to 
sustain a remarkably high degree of supersonic turbulence for a long time interval in the cloud. 

Finally, we would like to mention that soliton solutions may lead to a very different picture 
of the interstellar medium. The temperature fluctuations related to their wavy nature can create 
time-dependent patterns which can propagate and survive for long times due to their intrinsic 
stability. This may have some relevance for the interfaces producing the highly ionized species 
(C IV, N V and O VI) detected either in the disk and in the halo by absorption line measures. The 
study of extragalactic environments where two-phase media are found, {e.g., protogalactic clouds 
in a hot intergalactic medium, emitting clouds in the Broad Line Region of AGNs, cooling flows) 
constitutes an appealing application of the theory developed here. However, still more physical 
insight can be achieved by a dedicated analysis in which the three dimensional properties of soliton 
solutions are considered. 
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The solutions found for static CC fronts are intimately connected with the Field length. 
To see this, we recall that eq. (3.14) for k{T) = has the form 



For kt{T2) 7^ and Ct{T2) 7^ 0, this results in a finite size of the region with T ~ T2, of the order 
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and 




The fact that Ax(T2) comes out to be of the order of the Field length ip is not surprising. This is in 
agreement with the discussed interpretation of ip as the length over which cooling and conductive 
energy fluxes are able to establish equilibrium. 

The previous statement can be analogously demonstrated also in the case k{T) 7^ 0, in 
which conductivity provides heat transport over the whole system. At 77 = 770 (center), eq.(3.16) 
has the following solution 

q'^ = Qo-\ ^rio\iv-mf, 

where Qq = 2 C{ri)dri. It follows that 

i] = i]o±ql£%sm—, 
with = ^pivo)- Similarly, for rj = r]2 (saddle), the same equation has the solution 

with ip^ = ^F{ri2). 



APPENDIX B 

We show here that the flow described by eq. (3.39 is subsonic at ^ = +00. First we recall 

that 

|f/2| = |?7o|^^2^^|A.|; 

7 

If \Xs\ ~ A, then IC/2I ^ C2 and the flow is subsonic. This can be shown the following way. We can 
express the non-dimensional cooling function A in the equivalent form 

IpCmC £p ~ 

A —— A 

IclpiCm 4 

where Cm 7^ is the absolute value of the cooling function at a certain temperature Ti <T < T2; 
for example, we can assume that Cm is the maximum of \C\ in the interval (Ti, T2); ic = CiTc is the 

wavelength whose inverse frequence is equal to the cooling time; A = Cj Cm is the cooling function 
normalized to its maximal value Cm- The integral JJ^^ \(T)dC ~ A; this means that |As| '-^ ip/^c- 
Taking into account that k oc VT/{crein), where vt ~ ci is the thermal velocity, a^i is the cross 
section for clastic scattering of particles providing heat transfer, and C oc [a pinvj-) ^Epi^, where 
/S.E is tlicmean energy radiated by two colliding particles in inelastic processes, and is the mean 
probability for inelastic processes, we obtain Ip/lc ^ {Pin^E / kTi)^/^ . This estimate demonstrates 
that ip/ic is much smaller than unity because inelastic processes are sufficiently less frequent than 
elastic ones. For example, in a free-free radiating gas at T ~ 10^ K, pin^E/kT ^ vt/c where 
c is the light speed (Landau &: Lifshits 1975). These arguments demonstrate that evaporation or 
condensation in a stationary flow are always subsonic. 



REFERENCES 



Adams, F. C. & Fatuzzo, M. 1992, Univ. of Michigan, preprint 

Balbus, S. A. 1985, ApJ, 291, 518 

Balbus, S. A. 1986, ApJ, 381, 137 

Balbus, S. A. & McKee, C. F. 1982, ApJ, 252, 529 

Ballet, J., Arnaud, M. & Rothenfiug, R. 1986, A&A, 161, 12 

Burrows, D. N. & Mendenhall, J. A. 1991, Nature, 351, 629 

Bcgclman, M. C. & McKcc, C. F. 1990, ApJ, 358, 375 

Benjamin, A. 1972, Proc. Roy. Soc, A328, 153 

Bertoldi, F. 1989, ApJ, 346, 735 

Bertoldi, F. & McKee, C. F. 1990, ApJ, 354, 529 

Black, J. H. 1981, MNRAS, 197, 555 

Borkowski, K.J., Balbus, S. A., & Fristrom, C.C. 1990, ApJ, 355, 501 

Cowic L.L & McKee, C. F. 1977a, ApJ, 211, 135 

Dalgarno, A. & McCray, R. A. 1972, ARA&A, 10, 375 

Doroshkevich, A. G. & Zel'dovich, Ya. B. 1981, Sov. Phys. JETP, 53, 405 

Draine, B. T. & Giuliani, J. L. 1984, ApJ, 281, 690 

Elphick, C., Regev, O. & Spiegel, E. A. 1991, MNRAS, 250, 617 

Elphick, C., Regev, O. & Shaviv N. 1992, ApJ, 392, 106 

Ferrara, A. & Einaudi, G. 1992, ApJ, 395, 475 

Ferrara, A. 1993, ApJ, 407, 157 

Field, G. B. 1965, ApJ, 142, 531 

Field, G. B., Goldsmith D. W. & Habing, H. J. 1969, ApJ, 155, L149 
Giuliani, J. L. 1984, ApJ, 277, 605 

Guckenheimer, J. & Holmes, P. 1983, Nonlinear Oscillations, Dynamical Systems and Bifurcations 
of Vector Fields, (Berlin:Springer-Verlag) 

Hayashi, C. 1985, Nonlinear Oscillations in Physical Systems, (Princeton: Univ. Press) 

Houck, J. C. & Bregman, J. N. 1990, ApJ, 352, 506 

Kovalenko, I. 1993, J. Comp. Phys., submitted 

Kovalenko, I. & Shchekinov, Yu. 1992, Astron. Astophyis. Trans., 1, 129 

Landau, L.D. & Lifshits, E.M. 1975, The Classical Theory of Fields, ( Oxford :Pergamon) 

Li, F. & Ikeuchi, S. 1992, ApJ, 390, 405 

Lockman, F. J., Hobbs, L. M. & Shull, J. M. 1986, ApJ, 301, 380 
Lockman, F. J., & Gehman, C. S. 1991, ApJ, 382, 182 
McKee, C. F. & Cowie, L.L. 1977, ApJ, 215, 213 
McKee, C. F. & Ostriker, J. P. 1977, ApJ, 218, 148 
McKee, C. F. & Begelman, M. C. 1990, ApJ, 358, 392 
Penston, M. V. k Brown, F. E. 1970, MNRAS, 150, 373 



Rosner, R. & Tucker, W. H. 1989, ApJ 338, 761 

Raymond, J. C, Cox, D. P., & Smith, B. W. 1976, ApJ, 204, 290 

Sedov, L. I. 1957, Similarity and Dimensional Methods in Mechanics, (New York: Academic) 

Slavin, J. D., 1989, ApJ, 346, 718 

Sembach, K. R. & Savage, B. D. 1992, ApJSS, 83, 147 

Snowden, S. L., Mebold, U., Hirth, W., Herbstmeier, & Schmitt, J. H. M. M., 1991, Science, 252, 
1529 

Spitzer, L. 1962, Physics of Fully Ionized Gases, (New York:Interscience) 
Spitzer, L. 1990, ARA&A, 28, 71 
Tribble, P. C. 1989, MNRAS, 238, 1247 

Zel'dovich, Ya. B. & Kompaneets, A. S. 1955, Theory of Detonation, Moscow 

Zel'dovich, Ya. B. Sz Raizer 1965, Physics of Shock Waves and High Temperature Hydrodynamic 
Phenomena, (London: Academic Press) 

Zel'dovich, Ya. B. & Pikelner, S.B. 1969, Sov. Phys. JETP, 29, 170 



FIGURE CAPTIONS 



Figure 1 Different possible regimes for an inhomogeneous medium as a function of tlic ratio between 
the characteristic conductive and cooling times. The linear curve corresponds to eq. (2.13), the 
cubic curve corresponds to eq. (2.14), and the negative slope curve corresponds to eq.(2.14a); the 
vertical line corresponds to a fixed temperature of the hot phase; dashed regions correspond to 
states with large differences in the values of the characteristic times. 

Figure 2 Evolutionary stages of a CC front as a function of radius in units of i for the case if = oo. 
The parameters adopted are n = 1.0 cm~^, T = 10^ K, I = 3.68 pc, T = 3 x 10^ K. In each panel 
the flow Mach number {solid line), the density (dotted), and the log of temperature (dashed) and 
pressure (long- dashed), all normalized to the cloud values, are reported. Note that in the last four 
panels the density is divided by a factor of 10, as indicated by the label n/10. 

Figure 3 Evolution of the local saturation parameter cjt as a function of the cloud radius for the 
same case shown in Fig. 2. 

Figure 4 The same as Fig. 2, but (p = 2. 

Figure 5 The same as Fig. 2, density isocontours. 

Figure 6 Evolution of an isobaric thermal conduction front (numbers correspond to time elapsed 
from t = 0) with a real cooling function. Solid curves full solution, dashed curves no cooling 
included. 

Figure 7 Segments of phase trajectories near the stationary point ryi for a) > 0, b) Cr^ = 0, 
and c) jCrii = +oo. 

Figure 8 a) Phase diagram corresponding to the system (3.12), (3.13) ; b) eigensolution of (3.12), 
(3.13). 

Figure 9 Phase trajectories and eigenfunctions of eq. (3.12) with a) r(— oo) = T(+oo) = T2 (cold 
cloud immersed in a hot intercloud gas); b) r(— 00) = r(+oo) = Ti hot bubble in a cold infinite 
gas; c) periodic structure of cloud - intercloud gas. 

Figure 10 Van der Waals-type equation of state for interstellar gas. The dotted lines correspond 
to the intermediate region — 00 < ^ < +00. 

Figure 11a) Schematic dependence of the cooling function on ^: for an evaporation wave cooling 

function is asymmetric and heating dominates (Xg < 0); b) structure of the integral curves near 
the critical points; c) phase diagram for an evaporation CC front and d) for a condensation one. 

Figure 12 Asymmetric "potential" curve. The straight line shows the increase in potential energy 
due to a negative "friction" . 

Figure 13 Special case of a solitary a) condensation-evaporation and b) evaporation-condensation 
wave with fixed conductive heat fiux at the moving boundary S. 

Figure 14 a) Phase diagram and b) soliton-like profile of T for the CC front described by eq. (4.5) 
with 7 + 1/7 — 1 = F(T2); the thermally stable phase is located at T = r(iboo) = Ti. 



